Introduction

Depression is one of the most common mental health conditions worldwide. It can have significant consequences for individuals, communities, and wider society, including inmpacts on wellbeing and employment. One form of depression known to worsen during the winter months is Seasonal Affective Disorder (SAD) (springer/references).

The Covid-19 pandemic created major distruptions to daily life, including prolonged periods spent indoors due to lockdowns. Reduced exposure to natural light is a known risk factor for SAD(NHS), and the combination of winter daylight shortages and enforced indoor isolation may have contributed to increased levels of depressive symptoms during the pandemic. These changes may be reflected in pattern of antidepressant prescribing, particularly during winter periods.

The aim of this study is to investigate wether the COVID-19 pandemic influenced antidepressant prescribing and whether there are inequalities or geographic differences associated with this trend.

Methods

All data used in this study was obtained from publicly available NHS Scotland and Scottish Government sources.Antidepressant prescribing data were downloaded from the Data by Prescriber Location dataset, focusing specifically on winter periods (December to Februray) for the years 2019-2023, from the Public Health Scotland OpenData portal. Winter population estimates were taken from the National Records of Scotland census release.Socioeconomic information was sourced from the Scottish Index of Multiple Deprivation (SIMD 2020v2) datasets published by the Scottish Government. Spatial analysis was supported by incorporating NHS Scotland Health Board boundary shapefiles, which are publicly available through PHS and used to map regional variation in prescribing.

Antidepressant medications were identified using the BNF item code prefix “0403”, and all prescriptions with codes starting with “0403” were retained for analysis. These filtered data were then processed and summarised to generate visualisations of prescribing patterns.

For the purpose of this study, the COVID-19 period is defined as March 2020 to June 2021, corresponding to the main period of national lockdowns. This places Winter 2019 as pre-COVID, Winter 2020/2021 as during COVID, and Winter 2022 as post-COVID.

RESULTS

##Loading Necessary Packages
library(tidyverse) # <2>
library(here) # directory stucture
library(gt) # tables
library(janitor) # cleaning data
library(ggplot2) # plotting graph
library(sf) # to read in map data 
library(readxl) # to read in map data
library(plotly) # to make interactive
library(viridis) # colour palette for map
library(ggiraph)

2.tidyverse

# Loading multiple data files (18 CSVs) containing counts of medications dispensed in the community
files <- list.files(here("data", "winter_data"), pattern = "csv")
winter_data <- files %>% # Read all downloaded winter prescribing files stored in the 'winter_data' folder
  map_dfr(~read_csv(here("data", "winter_data", .))) %>% 
clean_names()#combine into one dataframe and clean column names
filtered_winter_data <- winter_data %>% 
filter(str_starts(bnf_item_code,"0403")) %>%  #Filter for antidepressants only (BNF item code starts with "0403" 
  mutate(year = as.numeric(substr(paid_date_month,1,4)), month = as.numeric(substr(paid_date_month,5,6))) %>% #Extract year and month from 'paid_date_month' for grouping into winter periods
  mutate(winter_year=case_when(month == 12 ~ year + 1, 
month %in% c(1,2) ~ year) )  #Create 'winter_year' column: December counts belong to the next year, Jan-Feb to current year
filtered_winter_data <- filtered_winter_data %>% 
  unite("healthboards",hbt2014,hbt,sep = "_") #Merge Health Board columns into a single column since they contain the same information
filtered_winter_data$healthboards <- gsub("[NA]","",filtered_winter_data$healthboards) 
    filtered_winter_data$healthboards <-
      gsub("_","",filtered_winter_data$healthboards)#Remove any remaining 'NA' strings and underscores to standardize Health Board codes

Figure 1

 #summarise total antidepressant prescriptions for each winter year
winter_years_data <- filtered_winter_data %>% 
  group_by(winter_year) %>% 
  summarise(total_items=sum(number_of_paid_items,na.rm = TRUE))
#plot line graph to visualise overall trend in total prescriptions across winter years
plot <- suppressWarnings( ggplot(winter_years_data,aes(x=winter_year,y=total_items)) + #line connecting total prescriptions per year
  geom_line(linewidth=0.7,colour = "blue") +
  geom_point(aes(text = paste0(
      "<b>Winter Year:</b> ",winter_year, "<br>",
      "<b>Total Items:</b> ",scales::comma(total_items)),
size = 1))+ #points with interactive tooltip showing year and total
  geom_vline(xintercept = 2020, 
linetype = "solid", colour = "brown") + geom_vline(xintercept = 2022, 
linetype = "solid", colour = "brown") + #reference lines to mark COVID-related winters
  scale_x_continuous(breaks=2017:2023) + #ensures all years appear on x-axis
  labs(title="Total Antidepressant Prescriptions During Winter Season",x="Winter Year",y="Total Antidepressant Prescriptions") +
  theme_minimal() + 
 theme(
    panel.background = element_rect(fill = "lightblue"),   # plot area colour
    plot.background = element_rect(fill = "white"), # whole plot colour
    panel.grid.major = element_line(colour="white",linewidth=0.5), # grid line colour and width
    panel.grid.minor = element_line(colour="white",linewidth=0.3) #gridline colour and width
  ))
  
ggplotly(plot,tooltip="text") #convert ggplot to interactive plotly object with tooltips

The line graph shows a steady upward trend in winter antidepressant prescribing across all years. Although the COVID-19 period (2020–2022) is highlighted, prescribing patterns do not markedly diverge from pre-pandemic trends, with increases during COVID appearing slightly less steep than earlier years. The y-axis is truncated to improve visibility of year-to-year changes. 

#Prepare Population
population <- readxl::read_excel(here("data","population.xlsx"), skip=10) %>% 
  clean_names() %>% 
  select(x2,all_people) %>% 
  filter(!is.na(all_people)) %>% 
  rename(h_bname = "x2",hb_population = "all_people") %>% 
filter(!str_detect(hb_population,"Cells"))

filtered3_winter_data <- filtered_winter_data %>% 
  group_by(healthboards,winter_year,gp_practice) %>% 
  summarise(paid_quantity = sum(paid_quantity,na.rm=TRUE)) 

SIMD <- readxl::read_excel(here("data","SIMD.xlsx")) %>% 
  clean_names() 

combined_pop_simd <- SIMD %>% 
  full_join(filtered3_winter_data,join_by(h_bcode == healthboards)) %>% 
  left_join(population,join_by(h_bname == h_bname))

Figure 2

#sum antidepressant prescriptions per health board and winter year (from 2019 onwards)
hb_prescribing <- filtered_winter_data %>%
  filter(winter_year >= 2019) %>%
  group_by(healthboards, winter_year) %>%
  summarise(total_paid = sum(paid_quantity, na.rm = TRUE), .groups="drop")
#clean SIMD data (keep unique health board codes and names)
simd_clean <- SIMD %>% select(h_bcode, h_bname) %>% distinct()
#combine with prescribing data with SIMD and population data,calculate prescriptions per head
combined <- hb_prescribing %>%
  left_join(simd_clean, join_by(healthboards == h_bcode)) %>%
  left_join(population, by = "h_bname") %>%
  mutate(quantity_per_head = total_paid / hb_population) %>%
  filter(!is.na(h_bname))
#read and clean NHS health board shapefile
NHS_healthboards <- st_read(here("data","Week6_NHS_HealthBoards_2019.shp")) %>%
  clean_names() %>% rename(h_bname = hb_name)
## Reading layer `Week6_NHS_healthboards_2019' from data source 
##   `/Users/olufimihanfaturoti/Year 3 Medicine/data-science/B251495/data/Week6_NHS_healthboards_2019.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 7564.996 ymin: 530635.8 xmax: 468754.8 ymax: 1218625
## Projected CRS: OSGB36 / British National Grid
#join spatial data with prescribing data, fill missing value with 0
retry_mapped_data <- NHS_healthboards %>%
  left_join(combined, by = "h_bname") %>%
  st_as_sf() %>%
  mutate(quantity_per_head = replace_na(quantity_per_head, 0))
#plot interactive map of antidepressant prescriptions per head
plot_map <- retry_mapped_data %>% 
  ggplot() +
  geom_sf_interactive(aes(fill = quantity_per_head,
        tooltip = paste0(h_bname,
                         "\nWinter Year: ", winter_year,
                         "\nItems per Head: ", round(quantity_per_head,2))),
    colour = "white",
    size = 0.1) +
  scale_fill_viridis_c(option = "mako", direction = -1, name = "Items per Head") + #blueish colour palette
  labs(title="Antidepressant Prescriptions per Head",
       subtitle="By Health Board and Winter Year") +
  facet_wrap(~winter_year) + # separate maps per year
  theme_void() +
  theme(strip.text = element_text(size=12, face="bold"),
        plot.title = element_text(face="bold", size=16),
        plot.subtitle = element_text(size=10))
#convert to interactive map
interactive_map <- girafe(ggobj = plot_map);interactive_map

results

Following from the overall increase in aantidepressant prescriptions shown in Figure 1, this map shows the geographic distribution of prescriotions and highlights both consistency and change. While the general north-south gradient in prescribing rates remained stable, the overall intensification post-2020 suggests a system wide response to increased levels depression. Interactive exploration of the map shows that health bards such as Lanarkshire and Glasgow experienced the largest increase in antidepressant prescriptions between the pre-COVID period(2019). The progression overall from 2020-2022 shows a clear deepinging of intesisty suggesting that regions with historically higher prescription rates saw even greater increases during and after the pandemic. This may reflect widening regional disparities. Overall, central Scotland, particularly the more densely populated urban belt, shows higher antidepressant prescriptions per head compared to the northern regions, suggesting that population density and urbanization may be associated with higher prescribing rates.

Figure 3

#read and clean winter population data
winter_population <- list.files(here("data", "winter_population"), pattern = "csv", full.names = TRUE) %>%
  map_dfr(read_csv) %>%
  clean_names() %>%
  select(date, practice_code, hb, all_ages) %>%
  mutate(winter_year = as.numeric(substr(date, 1, 4)),
    winter_year = case_when(
      winter_year == 2020 ~ 2019,
      winter_year == 2023 ~ 2022,
      TRUE ~ winter_year ),
    healthboards = hb) %>%
  select(-hb)
#summarise and filter combined population data for selected years
summary_winter_with_simd <- combined_pop_simd %>%
  filter(!is.na(simd2020v2_quintile), winter_year %in% c("2019","2022")) %>%
  group_by(simd2020v2_quintile, h_bname, gp_practice, winter_year) %>%
  summarise(avg_paid_over_winters = mean(paid_quantity), .groups = "drop") %>%
  rename(practice_code = gp_practice) %>%
  mutate(winter_year = as.numeric(winter_year))
# Join datasets
joined_data <- summary_winter_with_simd %>%
  full_join(winter_population, by = c("winter_year", "practice_code"))
# Calculate prescriptions per 1000
final_boxplot <- joined_data %>%
  group_by(simd2020v2_quintile, h_bname, practice_code, winter_year) %>%
  reframe(
    avg_per_1000 = (avg_paid_over_winters / all_ages) * 1000,
    year_label = factor(winter_year,
                        levels = c(2019, 2022),
                        labels = c("Pre-COVID (2019)", "Post-COVID (2022)"))
  )

#create boxplot of average prescriptions per 1000 by SIMD quintile, faceted by pre/post COVID
box2 <- final_boxplot %>%
  ggplot(aes(x = simd2020v2_quintile, y = avg_per_1000, fill = factor(simd2020v2_quintile))) +
  geom_boxplot() +
  scale_fill_viridis(discrete = TRUE, alpha = 0.6) +
  scale_y_continuous(labels = scales::label_comma()) +
  coord_cartesian(ylim = c(0, 70000))  + #scaled y axis graph to be able to visualise better
  theme_minimal() +
  theme(legend.position = "none") +
  facet_wrap(~year_label) +
  labs(
    title = "Average Winter Antidepressant Prescriptions Per by SIMD Quintile",
    x = "SIMD Quintile (1 = Most Deprived)",
    y = "Avg prescriptions per 1000"
  )

ggplotly(box2)

results

Contrary to expectations, the boxplot shows no strong socioeconomic gradient in antidepressant prescribing rates across SIMD quintiles. Although interactive inspection reveals slight differences in median values, the overall pattern remains relatively flat. Median prescribing rates roughly doubled across all SIMD groups between 2019 and 2022 (rising to around 28,000 per 1,000), but this shift appears uniform, suggesting a system-wide response to the pandemic’s mental-health impact rather than a disproportionate effect on specific socioeconomic groups. This contrasts with the clearer geographical differences observed in Figure 2. A flat pattern across deprivation levels could indicate equitable access to mental-health treatment. Alternatively, it may reflect potential under-prescribing in the most deprived areas (SIMD Quintile 1), where need would typically be expected to be higher. Both time periods show substantial upper-range outliers (50,000–70,000+ per 1,000), meaning some health boards or GP practices prescribe at rates two to three times above the median regardless of deprivation level or year. These persistent outliers likely point to local prescribing cultures, service availability, or population health differences, and warrant further investigation. While COVID-19 did not widen median socioeconomic inequalities in prescribing, it did increase variability within each SIMD group. The wider interquartile ranges and more extreme outliers in 2022 suggest growing disparities in how individual health boards or practices responded to mental-health demand during and after the pandemic, even when deprivation levels were similar.

Figure 4

# Create health board deprivation profile
hb_deprivation <- SIMD %>%
  group_by(h_bname) %>%
  summarise(
    # Population-weighted average quintile
    avg_quintile = weighted.mean(simd2020v2_quintile, w = population, na.rm = TRUE),
    
    # Percentage in most deprived quintiles (Q1 & Q2)
    pct_deprived = sum(population[simd2020v2_quintile <= 2], na.rm = TRUE) / 
                   sum(population, na.rm = TRUE) * 100,
    
    # Total population for reference
    total_pop = sum(population, na.rm = TRUE)
  ) %>%
  arrange(avg_quintile)  # Lower = more deprived

# Then join this to your prescribing data
final_lollipop_data <- combined %>% 
  filter(winter_year %in% c("2019","2022")) %>% 
  mutate(period = ifelse(winter_year == 2019, "pre", "post")) %>%
  left_join(hb_deprivation, by = "h_bname")

# Now create your lollipop plot ordered by deprivation
lollipop_plot <- final_lollipop_data %>% 
  group_by(h_bname, period) %>% 
  summarise(
    avg_prescribe = mean(total_paid),
    avg_quintile = first(avg_quintile),  # Take the pre-calculated value
    .groups = "drop"
  ) %>% 
  pivot_wider(names_from = period, values_from = avg_prescribe) %>%
  mutate(pct_change = ((post - pre) / pre) * 100) %>%
  arrange(avg_quintile) %>%  # Order by deprivation
  mutate(h_bname = factor(h_bname, levels = unique(h_bname)))

pct_lollipop <-lollipop_plot%>% 
  ggplot () +
  geom_segment(aes(x=h_bname, xend=h_bname, y=0, yend = pct_change), color = "skyblue")+
  geom_point(aes(x = h_bname, y = pct_change),color = "blue", size=3, alpha = 0.6) +
  theme_light() +
  coord_flip() +  
  theme(panel.grid.major.y = element_blank(),
    panel.border = element_blank(),
    axis.ticks.y = element_blank()) +
  labs( title = "Percentage Change in Prescribing per Health Board",
    x = "Health Board",
    y = "Percentage Change (%)")
   ggplotly(pct_lollipop)

Conclusion Limitations Next Steps